package triangleindex;
/*
 * @ this class is for dymaxion projection
 */

public class DymaxionProjection {

	/************************************************************************/
	/* NOTE: in C, array indexing starts with element zero(0). I choose */
	/* to start my array indexing with elemennt one(1) so all arrays */
	/* are defined one element longer than they need to be. */
	/************************************************************************/

	/************************************************************************/
	/* global variables accessible to all procedures */
	/************************************************************************/

	static float SQRT_3 = (float) Math.sqrt(3);
	static float SQRT_5 = (float) Math.sqrt(5);
	static float SQRT_8 = (float) Math.sqrt(8);
	static float SQRT_10 = (float) Math.sqrt(10);
	static float SQRT_15 = (float) Math.sqrt(15);

	double[] v_x = new double[13];
	double[] v_y = new double[13];
	double[] v_z = new double[13];
	double[] center_x = new double[21];
	double[] center_y = new double[21];
	double[] center_z = new double[21];
	double garc;
	double gt;
	double gdve;
	double gel;

	/****************************************/
	/* function definitions */
	/****************************************/
	double last_lng = 100000;
	double last_lat = 100000;
	point2d last_point; // ????

	public DymaxionProjection() {
		init_stuff();
	}

	public point2d convert_s_t_p_cache(double lng, double lat) {
		if (last_lng == lng && last_lat == lat) {
			return last_point;
		} else {
			last_lng = lng;
			last_lat = lat;
			last_point = convert_s_t_p(lng, lat);
			return last_point;
		}
	}

	public point2d convert_s_t_p(double lng, double lat) {
		/***********************************************************/
		/* This is the main control procedure. */
		/***********************************************************/

		/* Convert the given(long.,lat.) coordinate into spherical */
		/* polar coordinates(r, theta, phi) with radius=1. */
		/* Angles are given in radians, NOT degrees. */

		double[] sc = conv_ll_t_sc(lng, lat);

		/* convert the spherical polar coordinates into cartesian */
		/* (x, y, z) coordinates. */

		point3d h = s_to_c(sc[0], sc[1]); // sc[0] = theta, sc[1] = phi

		/* determine which of the 20 spherical icosahedron triangles */
		/* the given point is in and the LCD triangle. */

		int[] info = s_tri_info(h.x, h.y, h.z);

	//	System.out.println("info[0]" + info[0] + "info[1]" + info[1]);
		/* Determine the corresponding Fuller map plane(x, y) point */

		return dymax_point(info[0], info[1], h.x, h.y, h.z);
		
	}

	public point2d convert_i_t_p(int i, int[] t) {
		double x = 0;
		double y = 0;
		double z = 0;

		// We need to 'nudge' the point a little bit into the triangle
		// Hence we do a weighted average with point i having a massive weight
		for (int j = 0; j < 3; j++) {
			if (t[j] == i) {
				x += v_x[i] * 0.9999;
				y += v_y[i] * 0.9999;
				z += v_z[i] * 0.9999;
			} else {
				x += v_x[t[j]] * 0.00005;
				y += v_y[t[j]] * 0.00005;
				z += v_z[t[j]] * 0.00005;
			}
		}

		int[] info = s_tri_info(x, y, z);
		/* Determine the corresponding Fuller map plane(x, y) point */

		return dymax_point(info[0], info[1], x, y, z);
	}

	public double[] conv_ll_t_sc(double lng, double lat) {
		/* convert(long., lat.) point into spherical polar coordinates */
		/* with r=radius=1. Angles are given in radians. */
		double[] sc = new double[2];
		double h_theta, h_phi;

		h_theta = 90.0 - lat;
		h_phi = lng;
		if (lng < 0.0) {
			h_phi = lng + 360.0;
		}
		// sc[0] = theta sc[1]=phi
		sc[0] = radians(h_theta);
		sc[1] = radians(h_phi);

		return sc;
	} /* end conv_ll_t_sc */

	// tranform degree to PI
	public double radians(double degrees) {
		/* convert angles in degrees into angles in radians */
		return (Math.PI * degrees / 180);
	} /* end of radians function */

	public void init_stuff() {
		/* initializes the global variables which includes the */
		/* vertix coordinates and mid-face coordinates. */

		double i, hold_x, hold_y, hold_z, magn;
		double theta, phi;

		/* Cartesian coordinates for the 12 vertices of icosahedron */

		v_x[1] = 0.420152426708710003;
		v_y[1] = 0.078145249402782959;
		v_z[1] = 0.904082550615019298;

		v_x[2] = 0.995009439436241649;
		v_y[2] = -0.091347795276427931;
		v_z[2] = 0.040147175877166645;

		v_x[3] = 0.518836730327364437;
		v_y[3] = 0.835420380378235850;
		v_z[3] = 0.181331837557262454;

		v_x[4] = -0.414682225320335218;
		v_y[4] = 0.655962405434800777;
		v_z[4] = 0.630675807891475371;

		v_x[5] = -0.515455959944041808;
		v_y[5] = -0.381716898287133011;
		v_z[5] = 0.767200992517747538;

		v_x[6] = 0.355781402532944713;
		v_y[6] = -0.843580002466178147;
		v_z[6] = 0.402234226602925571;

		v_x[7] = 0.414682225320335218;
		v_y[7] = -0.655962405434800777;
		v_z[7] = -0.630675807891475371;

		v_x[8] = 0.515455959944041808;
		v_y[8] = 0.381716898287133011;
		v_z[8] = -0.767200992517747538;

		v_x[9] = -0.355781402532944713;
		v_y[9] = 0.843580002466178147;
		v_z[9] = -0.402234226602925571;

		v_x[10] = -0.995009439436241649;
		v_y[10] = 0.091347795276427931;
		v_z[10] = -0.040147175877166645;

		v_x[11] = -0.518836730327364437;
		v_y[11] = -0.835420380378235850;
		v_z[11] = -0.181331837557262454;

		v_x[12] = -0.420152426708710003;
		v_y[12] = -0.078145249402782959;
		v_z[12] = -0.904082550615019298;

		/* now calculate mid face coordinates */

		hold_x = (v_x[1] + v_x[2] + v_x[3]) / 3.0;
		hold_y = (v_y[1] + v_y[2] + v_y[3]) / 3.0;
		hold_z = (v_z[1] + v_z[2] + v_z[3]) / 3.0;
		magn = Math.sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
		center_x[1] = hold_x / magn;
		center_y[1] = hold_y / magn;
		center_z[1] = hold_z / magn;

		hold_x = (v_x[1] + v_x[3] + v_x[4]) / 3.0;
		hold_y = (v_y[1] + v_y[3] + v_y[4]) / 3.0;
		hold_z = (v_z[1] + v_z[3] + v_z[4]) / 3.0;
		magn = Math.sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
		center_x[2] = hold_x / magn;
		center_y[2] = hold_y / magn;
		center_z[2] = hold_z / magn;

		hold_x = (v_x[1] + v_x[4] + v_x[5]) / 3.0;
		hold_y = (v_y[1] + v_y[4] + v_y[5]) / 3.0;
		hold_z = (v_z[1] + v_z[4] + v_z[5]) / 3.0;
		magn = Math.sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
		center_x[3] = hold_x / magn;
		center_y[3] = hold_y / magn;
		center_z[3] = hold_z / magn;

		hold_x = (v_x[1] + v_x[5] + v_x[6]) / 3.0;
		hold_y = (v_y[1] + v_y[5] + v_y[6]) / 3.0;
		hold_z = (v_z[1] + v_z[5] + v_z[6]) / 3.0;
		magn = Math.sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
		center_x[4] = hold_x / magn;
		center_y[4] = hold_y / magn;
		center_z[4] = hold_z / magn;

		hold_x = (v_x[1] + v_x[2] + v_x[6]) / 3.0;
		hold_y = (v_y[1] + v_y[2] + v_y[6]) / 3.0;
		hold_z = (v_z[1] + v_z[2] + v_z[6]) / 3.0;
		magn = Math.sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
		center_x[5] = hold_x / magn;
		center_y[5] = hold_y / magn;
		center_z[5] = hold_z / magn;

		hold_x = (v_x[2] + v_x[3] + v_x[8]) / 3.0;
		hold_y = (v_y[2] + v_y[3] + v_y[8]) / 3.0;
		hold_z = (v_z[2] + v_z[3] + v_z[8]) / 3.0;
		magn = Math.sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
		center_x[6] = hold_x / magn;
		center_y[6] = hold_y / magn;
		center_z[6] = hold_z / magn;

		hold_x = (v_x[8] + v_x[3] + v_x[9]) / 3.0;
		hold_y = (v_y[8] + v_y[3] + v_y[9]) / 3.0;
		hold_z = (v_z[8] + v_z[3] + v_z[9]) / 3.0;
		magn = Math.sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
		center_x[7] = hold_x / magn;
		center_y[7] = hold_y / magn;
		center_z[7] = hold_z / magn;

		hold_x = (v_x[9] + v_x[3] + v_x[4]) / 3.0;
		hold_y = (v_y[9] + v_y[3] + v_y[4]) / 3.0;
		hold_z = (v_z[9] + v_z[3] + v_z[4]) / 3.0;
		magn = Math.sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
		center_x[8] = hold_x / magn;
		center_y[8] = hold_y / magn;
		center_z[8] = hold_z / magn;

		hold_x = (v_x[10] + v_x[9] + v_x[4]) / 3.0;
		hold_y = (v_y[10] + v_y[9] + v_y[4]) / 3.0;
		hold_z = (v_z[10] + v_z[9] + v_z[4]) / 3.0;
		magn = Math.sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
		center_x[9] = hold_x / magn;
		center_y[9] = hold_y / magn;
		center_z[9] = hold_z / magn;

		hold_x = (v_x[5] + v_x[10] + v_x[4]) / 3.0;
		hold_y = (v_y[5] + v_y[10] + v_y[4]) / 3.0;
		hold_z = (v_z[5] + v_z[10] + v_z[4]) / 3.0;
		magn = Math.sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
		center_x[10] = hold_x / magn;
		center_y[10] = hold_y / magn;
		center_z[10] = hold_z / magn;

		hold_x = (v_x[5] + v_x[11] + v_x[10]) / 3.0;
		hold_y = (v_y[5] + v_y[11] + v_y[10]) / 3.0;
		hold_z = (v_z[5] + v_z[11] + v_z[10]) / 3.0;
		magn = Math.sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
		center_x[11] = hold_x / magn;
		center_y[11] = hold_y / magn;
		center_z[11] = hold_z / magn;

		hold_x = (v_x[5] + v_x[6] + v_x[11]) / 3.0;
		hold_y = (v_y[5] + v_y[6] + v_y[11]) / 3.0;
		hold_z = (v_z[5] + v_z[6] + v_z[11]) / 3.0;
		magn = Math.sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
		center_x[12] = hold_x / magn;
		center_y[12] = hold_y / magn;
		center_z[12] = hold_z / magn;

		hold_x = (v_x[11] + v_x[6] + v_x[7]) / 3.0;
		hold_y = (v_y[11] + v_y[6] + v_y[7]) / 3.0;
		hold_z = (v_z[11] + v_z[6] + v_z[7]) / 3.0;
		magn = Math.sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
		center_x[13] = hold_x / magn;
		center_y[13] = hold_y / magn;
		center_z[13] = hold_z / magn;

		hold_x = (v_x[7] + v_x[6] + v_x[2]) / 3.0;
		hold_y = (v_y[7] + v_y[6] + v_y[2]) / 3.0;
		hold_z = (v_z[7] + v_z[6] + v_z[2]) / 3.0;
		magn = Math.sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
		center_x[14] = hold_x / magn;
		center_y[14] = hold_y / magn;
		center_z[14] = hold_z / magn;

		hold_x = (v_x[8] + v_x[7] + v_x[2]) / 3.0;
		hold_y = (v_y[8] + v_y[7] + v_y[2]) / 3.0;
		hold_z = (v_z[8] + v_z[7] + v_z[2]) / 3.0;
		magn = Math.sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
		center_x[15] = hold_x / magn;
		center_y[15] = hold_y / magn;
		center_z[15] = hold_z / magn;

		hold_x = (v_x[12] + v_x[9] + v_x[8]) / 3.0;
		hold_y = (v_y[12] + v_y[9] + v_y[8]) / 3.0;
		hold_z = (v_z[12] + v_z[9] + v_z[8]) / 3.0;
		magn = Math.sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
		center_x[16] = hold_x / magn;
		center_y[16] = hold_y / magn;
		center_z[16] = hold_z / magn;

		hold_x = (v_x[12] + v_x[9] + v_x[10]) / 3.0;
		hold_y = (v_y[12] + v_y[9] + v_y[10]) / 3.0;
		hold_z = (v_z[12] + v_z[9] + v_z[10]) / 3.0;
		magn = Math.sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
		center_x[17] = hold_x / magn;
		center_y[17] = hold_y / magn;
		center_z[17] = hold_z / magn;

		hold_x = (v_x[12] + v_x[11] + v_x[10]) / 3.0;
		hold_y = (v_y[12] + v_y[11] + v_y[10]) / 3.0;
		hold_z = (v_z[12] + v_z[11] + v_z[10]) / 3.0;
		magn = Math.sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
		center_x[18] = hold_x / magn;
		center_y[18] = hold_y / magn;
		center_z[18] = hold_z / magn;

		hold_x = (v_x[12] + v_x[11] + v_x[7]) / 3.0;
		hold_y = (v_y[12] + v_y[11] + v_y[7]) / 3.0;
		hold_z = (v_z[12] + v_z[11] + v_z[7]) / 3.0;
		magn = Math.sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
		center_x[19] = hold_x / magn;
		center_y[19] = hold_y / magn;
		center_z[19] = hold_z / magn;

		hold_x = (v_x[12] + v_x[8] + v_x[7]) / 3.0;
		hold_y = (v_y[12] + v_y[8] + v_y[7]) / 3.0;
		hold_z = (v_z[12] + v_z[8] + v_z[7]) / 3.0;
		magn = Math.sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
		center_x[20] = hold_x / magn;
		center_y[20] = hold_y / magn;
		center_z[20] = hold_z / magn;

		garc = 2.0 * Math.asin(Math.sqrt(5 - SQRT_5) / SQRT_10);
		gt = garc / 2.0;

		gdve = Math.sqrt(3 + SQRT_5) / Math.sqrt(5 + SQRT_5);

		gel = SQRT_8 / Math.sqrt(5 + SQRT_5);
	} /* end of int_stuff procedure */

	public point3d s_to_c(double theta, double phi) {
		/* Covert spherical polar coordinates to cartesian coordinates. */
		/* The angles are given in radians. */
		point3d c = new point3d();
		c.x = Math.sin(theta) * Math.cos(phi);
		c.y = Math.sin(theta) * Math.sin(phi);
		c.z = Math.cos(theta);

		return c;
	} /* end s_to_c */

	public point2d c_to_s(double x, double y, double z) {
		/* convert cartesian coordinates into spherical polar coordinates. */
		/* The angles are given in radians. */
		point2d s = new point2d();
		double a = 0;

		if (x > 0.0 && y > 0.0) {
			a = radians(0.0);
		}
		if (x < 0.0 && y > 0.0) {
			a = radians(180.0);
		}
		if (x < 0.0 && y < 0.0) {
			a = radians(180.0);
		}
		if (x > 0.0 && y < 0.0) {
			a = radians(360.0);
		}

		// s.x = s.lng; s.y =s.lat;
		s.y = Math.acos(z);
		if (x == 0.0 && y > 0.0) {
			s.x = radians(90.0);
		}
		if (x == 0.0 && y < 0.0) {
			s.x = radians(270.0);
		}
		if (x > 0.0 && y == 0.0) {
			s.x = radians(0.0);
		}
		if (x < 0.0 && y == 0.0) {
			s.x = radians(180.0);
		}
		if (x != 0.0 && y != 0.0) {
			s.x = Math.atan(y / x) + a;
		}

		return s;
	} /* end c_to_s */

	public int[] s_tri_info(double x, double y, double z) {
		/* Determine which triangle and LCD triangle the point is in. */

		double h_dist1, h_dist2, h_dist3, h1, h2, h3; // double
		int i, h_tri, h_lcd = 0; // int
		int v1 = 0, v2 = 0, v3 = 0; // int

		int[] info = new int[2]; // info[0] = tri info[1] = lcd

		h_tri = 0;
		h_dist1 = 9999.0;

		/* Which triangle face center is the closest to the given point */
		/* is the triangle in which the given point is in. */

		for (i = 1; i <= 20; i = i + 1) {
			h1 = center_x[i] - x;
			h2 = center_y[i] - y;
			h3 = center_z[i] - z;
			h_dist2 = Math.sqrt(h1 * h1 + h2 * h2 + h3 * h3);
			if (h_dist2 < h_dist1) {
				h_tri = i;
				h_dist1 = h_dist2;
			} /* end the if statement */
		} /* end the for statement */

		// info[0] = info.tri, indicate which triangle
		info[0] = h_tri;

		/* Now the LCD triangle is determined. */

		switch (h_tri) {
		case 1:
			v1 = 1;
			v2 = 3;
			v3 = 2;
			break;
		case 2:
			v1 = 1;
			v2 = 4;
			v3 = 3;
			break;
		case 4:
			v1 = 1;
			v2 = 5;
			v3 = 4;
			break;
		case 3:
			v1 = 1;
			v2 = 6;
			v3 = 5;
			break;
		case 5:
			v1 = 1;
			v2 = 2;
			v3 = 6;
			break;
		case 6:
			v1 = 2;
			v2 = 3;
			v3 = 8;
			break;
		case 7:
			v1 = 3;
			v2 = 9;
			v3 = 8;
			break;
		case 8:
			v1 = 3;
			v2 = 4;
			v3 = 9;
			break;
		case 9:
			v1 = 4;
			v2 = 10;
			v3 = 9;
			break;
		case 10:
			v1 = 4;
			v2 = 5;
			v3 = 10;
			break;
		case 11:
			v1 = 5;
			v2 = 11;
			v3 = 10;
			break;
		case 12:
			v1 = 5;
			v2 = 6;
			v3 = 11;
			break;
		case 13:
			v1 = 6;
			v2 = 7;
			v3 = 11;
			break;
		case 14:
			v1 = 2;
			v2 = 7;
			v3 = 6;
			break;
		case 15:
			v1 = 2;
			v2 = 8;
			v3 = 7;
			break;
		case 16:
			v1 = 8;
			v2 = 9;
			v3 = 12;
			break;
		case 17:
			v1 = 9;
			v2 = 10;
			v3 = 12;
			break;
		case 18:
			v1 = 10;
			v2 = 11;
			v3 = 12;
			break;
		case 19:
			v1 = 11;
			v2 = 7;
			v3 = 12;
			break;
		case 20:
			v1 = 8;
			v2 = 12;
			v3 = 7;
			break;
		} /* end of switch statement */

		h1 = x - v_x[v1];
		h2 = y - v_y[v1];
		h3 = z - v_z[v1];
		h_dist1 = Math.sqrt(h1 * h1 + h2 * h2 + h3 * h3);

		h1 = x - v_x[v2];
		h2 = y - v_y[v2];
		h3 = z - v_z[v2];
		h_dist2 = Math.sqrt(h1 * h1 + h2 * h2 + h3 * h3);

		h1 = x - v_x[v3];
		h2 = y - v_y[v3];
		h3 = z - v_z[v3];
		h_dist3 = Math.sqrt(h1 * h1 + h2 * h2 + h3 * h3);

		if ((h_dist1 <= h_dist2) && (h_dist2 <= h_dist3)) {
			h_lcd = 1;
		}
		if ((h_dist1 <= h_dist3) && (h_dist3 <= h_dist2)) {
			h_lcd = 6;
		}
		if ((h_dist2 <= h_dist1) && (h_dist1 <= h_dist3)) {
			h_lcd = 2;
		}
		if ((h_dist2 <= h_dist3) && (h_dist3 <= h_dist1)) {
			h_lcd = 3;
		}
		if ((h_dist3 <= h_dist1) && (h_dist1 <= h_dist2)) {
			h_lcd = 5;
		}
		if ((h_dist3 <= h_dist2) && (h_dist2 <= h_dist1)) {
			h_lcd = 4;
		}

		info[1] = h_lcd;

		return info;
	} /* end s_tri_info */

	public point2d dymax_point(int tri, int lcd, double x, double y, double z) {
		int axis, v1 = 0; // int
		point2d h;

		double gs; // double
		double gx, gy, gz, ga1, ga2, ga3, ga1p, ga2p, ga3p, gxp, gyp, gzp; // double

		/* In order to rotate the given point into the template spherical */
		/* triangle, we need the spherical polar coordinates of the center */
		/* of the face and one of the face vertices. So set up which vertex */
		/* to use. */

		switch (tri) {
		case 1:
			v1 = 1;
			break;
		case 2:
			v1 = 1;
			break;
		case 3:
			v1 = 1;
			break;
		case 4:
			v1 = 1;
			break;
		case 5:
			v1 = 1;
			break;
		case 6:
			v1 = 2;
			break;
		case 7:
			v1 = 3;
			break;
		case 8:
			v1 = 3;
			break;
		case 9:
			v1 = 4;
			break;
		case 10:
			v1 = 4;
			break;
		case 11:
			v1 = 5;
			break;
		case 12:
			v1 = 5;
			break;
		case 13:
			v1 = 6;
			break;
		case 14:
			v1 = 2;
			break;
		case 15:
			v1 = 2;
			break;
		case 16:
			v1 = 8;
			break;
		case 17:
			v1 = 9;
			break;
		case 18:
			v1 = 10;
			break;
		case 19:
			v1 = 11;
			break;
		case 20:
			v1 = 8;
			break;
		} /* end of switch statement */

		point3d h0 = new point3d();
		// h0[0] = h0.x; h0[1] = h0.y ; h0[2] = h0.z
		h0.x = x;
		h0.y = y;
		h0.z = z;

		point3d h1 = new point3d();
		h1.x = v_x[v1];
		h1.y = v_y[v1];
		h1.z = v_z[v1];

		h = c_to_s(center_x[tri], center_y[tri], center_z[tri]);

		axis = 3;
		rotate3d(axis, h.x, h0);
		rotate3d(axis, h.x, h1);

		axis = 2;
		rotate3d(axis, h.y, h0);
		rotate3d(axis, h.y, h1);

		h = c_to_s(h1.x, h1.y, h1.z);
		h.x = h.x - radians(90.0);

		axis = 3;
		rotate3d(axis, h.x, h0);

		/* exact transformation equations */

		gz = Math.sqrt(1 - h0.x * h0.x - h0.y * h0.y);
		gs = Math.sqrt(5 + 2 * SQRT_5) / (gz * SQRT_15);

		gxp = h0.x * gs;
		gyp = h0.y * gs;

		ga1p = 2.0 * gyp / SQRT_3 + (gel / 3.0);
		ga2p = gxp - (gyp / SQRT_3) + (gel / 3.0);
		ga3p = (gel / 3.0) - gxp - (gyp / SQRT_3);

		ga1 = gt + Math.atan((ga1p - 0.5 * gel) / gdve);
		ga2 = gt + Math.atan((ga2p - 0.5 * gel) / gdve);
		ga3 = gt + Math.atan((ga3p - 0.5 * gel) / gdve);

		gx = 0.5 * (ga2 - ga3);

		gy = (1.0 / (2.0 * SQRT_3)) * (2 * ga1 - ga2 - ga3);

		/* Re-scale so plane triangle edge length is 1. */
		point2d pt = new point2d();
		pt.x = gx / garc;
		pt.y = gy / garc;

		/* rotate and translate to correct position */
		point2d p2d = new point2d();

		switch (tri) {
		case 1:
			rotate2d(240.0, pt);
			p2d.x = pt.x + 2.0;
			p2d.y = pt.y + 7.0 / (2.0 * SQRT_3);
			break;
		case 2:
			rotate2d(300.0, pt);
			p2d.x = pt.x + 2.0;
			p2d.y = pt.y + 5.0 / (2.0 * SQRT_3);
			break;
		case 3:
			rotate2d(0.0, pt);
			p2d.x = pt.x + 2.5;
			p2d.y = pt.y + 2.0 / SQRT_3;
			break;
		case 4:
			rotate2d(60.0, pt);
			p2d.x = pt.x + 3.0;
			p2d.y = pt.y + 5.0 / (2.0 * SQRT_3);
			break;
		case 5:
			rotate2d(180.0, pt);
			p2d.x = pt.x + 2.5;
			p2d.y = pt.y + 4.0 * SQRT_3 / 3.0;
			break;
		case 6:
			rotate2d(300.0, pt);
			p2d.x = pt.x + 1.5;
			p2d.y = pt.y + 4.0 * SQRT_3 / 3.0;
			break;
		case 7:
			rotate2d(300.0, pt);
			p2d.x = pt.x + 1.0;
			p2d.y = pt.y + 5.0 / (2.0 * SQRT_3);
			break;
		case 8:
			rotate2d(0.0, pt);
			p2d.x = pt.x + 1.5;
			p2d.y = pt.y + 2.0 / SQRT_3;
			break;
		case 9:
			if (lcd > 2) {
				rotate2d(300.0, pt);
				p2d.x = pt.x + 1.5;
				p2d.y = pt.y + 1.0 / SQRT_3;
			} else {
				rotate2d(0.0, pt);
				p2d.x = pt.x + 2.0;
				p2d.y = pt.y + 1.0 / (2.0 * SQRT_3);
			}
			break;

		case 10:
			rotate2d(60.0, pt);
			p2d.x = pt.x + 2.5;
			p2d.y = pt.y + 1.0 / SQRT_3;
			break;
		case 11:
			rotate2d(60.0, pt);
			p2d.x = pt.x + 3.5;
			p2d.y = pt.y + 1.0 / SQRT_3;
			break;
		case 12:
			rotate2d(120.0, pt);
			p2d.x = pt.x + 3.5;
			p2d.y = pt.y + 2.0 / SQRT_3;
			break;
		case 13:
			rotate2d(60.0, pt);
			p2d.x = pt.x + 4.0;
			p2d.y = pt.y + 5.0 / (2.0 * SQRT_3);
			break;
		case 14:
			rotate2d(0.0, pt);
			p2d.x = pt.x + 4.0;
			p2d.y = pt.y + 7.0 / (2.0 * SQRT_3);
			break;
		case 15:
			rotate2d(0.0, pt);
			p2d.x = pt.x + 5.0;
			p2d.y = pt.y + 7.0 / (2.0 * SQRT_3);
			break;
		case 16:
			if (lcd < 4) {
				rotate2d(60.0, pt);
				p2d.x = pt.x + 0.5;
				p2d.y = pt.y + 1.0 / SQRT_3;
			} else {
				rotate2d(0.0, pt);
				p2d.x = pt.x + 5.5;
				p2d.y = pt.y + 2.0 / SQRT_3;
			}
			break;
		case 17:
			rotate2d(0.0, pt);
			p2d.x = pt.x + 1.0;
			p2d.y = pt.y + 1.0 / (2.0 * SQRT_3);
			break;
		case 18:
			rotate2d(120.0, pt);
			p2d.x = pt.x + 4.0;
			p2d.y = pt.y + 1.0 / (2.0 * SQRT_3);
			break;
		case 19:
			rotate2d(120.0, pt);
			p2d.x = pt.x + 4.5;
			p2d.y = pt.y + 2.0 / SQRT_3;
			break;
		case 20:
			rotate2d(300.0, pt);
			p2d.x = pt.x + 5.0;
			p2d.y = pt.y + 5.0 / (2.0 * SQRT_3);
			break;

		} /* end switch statement */

		return p2d;
	} /* end of dymax_point */

	public point2d rotate2d(double angle, point2d p2d) {
		/* Rotate the point to correct orientation in XY-plane. */

		double ha, hx, hy; // double

		ha = radians(angle);
		hx = p2d.x;
		hy = p2d.y;
		p2d.x = hx * Math.cos(ha) - hy * Math.sin(ha);
		p2d.y = hx * Math.sin(ha) + hy * Math.cos(ha);

		return p2d;
	} /* end rotate procedure */

	public point3d rotate3d(int axis, double alpha, point3d p) {
		/* Rotate a 3-D point about the specified axis. */

		double a = p.x;
		double b = p.y;
		double c = p.z;

		if (axis == 1) {
			p.y = b * Math.cos(alpha) + c * Math.sin(alpha);
			p.z = c * Math.cos(alpha) - b * Math.sin(alpha);
		}

		if (axis == 2) {
			p.x = a * Math.cos(alpha) - c * Math.sin(alpha);
			p.z = a * Math.sin(alpha) + c * Math.cos(alpha);
		}

		if (axis == 3) {
			p.x = a * Math.cos(alpha) + b * Math.sin(alpha);
			p.y = b * Math.cos(alpha) - a * Math.sin(alpha);
		}

		return p;
	} /* end of r2 */

}

class point3d {

	double x, y, z;

	public point3d() {
		this.x = 0;
		this.y = 0;
		this.z = 0;
	}

	public point3d(double x, double y, double z) {
		this.x = x;
		this.y = y;
		this.z = z;

	}

}

